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Abstract. We present numerical simulations of a self-sustaining magnetic field in a differentially rotating non-convective 
stellar interior. A weak initial field is wound up by the differential rotation; the resulting azimuthal field becomes unstable and 
produces a new meridional field component, which is then wound up anew, thus completing the 'dynamo loop' . This effect is 
observed both with and without a stable stratification. A self-sustained field is actually obtained more easily in the presence of 
a stable stratification. The results confirm the analytical expectations of the role of Tayler instability. 
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1. Introduction 

It has long been known that magnetic fields can be generated 
by a dynamo operating in the convective zone of a differentially 
rotating star (e.g. Parke r 19791 . A toroidal field is produced by 
winding-up of the poloidal (meridional) component, and the 
bubbles of gas moving upwards and downwards move perpen- 
dicular to these toroidal field lines, bending them and creating a 
new poloidal component, closing the 'dynamo loop'. This type 
of dynamo has been the subject of extensive work over several 
decades. A process like this has been held responsible for the 
magnetism of stars with convective envelopes like the Sun. 

Whether or not such direct causation by convection is ac- 
tually the correct explanation of dynamos like the solar cycle 1 , 
the predominance of this view has obscured the fact that self- 
sustained magnetic fields do not require the presence of convec- 
tion or other imposed small scale velocity fields. A magnetic 
field can produce small-scale perturbations from its own insta- 
bility, without recourse to externally imposed perturbations. A 
well-known example in the context of accretion discs, where 
a dynamo was produced when differential rotation wound up 
a field which was then subject to magnetohydrodynamic insta- 
bility (Hawley et al. 1996). 

The same principle was applied by Spruit (2002) to differ- 
entially rotating stars. In this scenario, a toroidal field is wound 
up by differential rotation from a weak seed field. The field re- 
mains predominantly toroidal, subject to decay by instability 
of the field, but is continuously regenerated by the winding- 
up of irregularities produced by the instability. This scenario 
has been applied in stellar evolution calculations of the inter- 
nal rotation of massive stars by Heger et al. (2003) and Maeder 
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1 It is in fact more likely that buoyant instability of the magnetic 
field, rather than convection, is the key process in the solar cycle, cf. 
the discussion in Spruit (1999). 



& Meynet (2003). The process is conceptually similar to the 
small-scale self-sustained fields found in MHD simulations of 
accretion discs, but operates on a different form of MHD in- 
stability (pinch-type or Tayler instability as opposed to magne- 
torotational, cf. Spruit 1999). 

A self-sustained field of this type could have important im- 
plications for not only the magnetism of a star, but also for 
the transfer of angular momentum. Differential rotation is cre- 
ated, when the star is formed, as a consequence of conserva- 
tion of angular momentum when parts of it contract or ex- 
pand, and through angular momentum loss through a stellar 
wind ('magnetic braking'). In the absence of a magnetic field, 
kinetic viscosity would eventually damp differential rotation, 
but only on a time-scale much longer than the lifetime of the 
star. If a weak magnetic field were present in a star with in- 
finite conductivity, the field would be wound up, its Lorentz 
force exerting a force back on the gas, tending to slow the dif- 
ferential rotation. If we assume that no magnetohydrodynamic 
instabilities were present, the energy of the field would be- 
come eventually comparable to the kinetic energy of the dif- 
ferential rotation, and the field would exert a force on the gas 
strong enough to reverse the differential rotation. Oscillations 
would follow, with energy continuously being transferred from 
kinetic to magnetic and back again (Meste l 1953ft . Finite con- 
ductivity would have the effect of damping these oscillations. 
However, if the magnetic field became unstable, as we expect 
it to, the energy of the field need never reach a level compara- 
ble to the kinetic energy and the direction of differential ro- 
tation would never be reversed. Instead, differential rotation 
would gradually be slowed, and the magnetic field held at a 
low steady-state level. This could be what has happened in the 
radiative core of the Sun, explaining the near-uniform rotation 
there JSchou et al. 19981 ICharbonneau et al. 19991 . 

While such a dynamo process is plausible, it has so far been 
described only in terms of an elementary scaling model (Spruit 
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2002). With the calculations presented here we verify, first of 
all, the existence of a self-sustained field generation process. 
At the next level, the goal is to compare the properties of the 
numerical results with the predictions of the model, in partic- 
ular concerning the central role of Tayler instability. The field 
strength resulting from the dynamo process depends on the bal- 
ance between the decay of the toroidal field and the winding 
up of irregularities. The analytic model only presents the basic 
scaling of this balance with parameters like the strength of the 
differential rotation; the actual values of the coefficients in this 
scaling have to be found by numerical simulations. 

1.1. Instability of a toroidal field 

Tayler (1973) and Acheson (1978) looked at toroidal fields in 
stars, that is, fields that have only an azimuthal component 
in a cylindrical coordinate frame (zu, <fi, z) with the origin at 
the centre of the star. With the energy method, Tayler derived 
necessary and sufficient stability conditions in adiabatic con- 
ditions (no viscosity, thermal diffusion or magnetic diffusion). 
The main conclusion was that such purely toroidal fields are 
always unstable at some place in the star, in particular to per- 
turbations of the m — 1 form, and that stability at any par- 
ticular place does not depend on field strength but only on the 
geometry of the field configuration. An important corollary of 
the results of Tayler (1973, esp. the Appendix) was the proof 
that instability is local in meridional planes. If the necessary 
and sufficient condition for instability is satisfied at any point 
(m, z), there is an unstable eigenfunction that will fit inside 
an infinitesimal environment of this point. The instability is al- 
ways global in the azimuthal direction, however. The instability 
takes place in the form of a low-azimuthal order displacement 
in a ring around the star. Connected with this is the fact that the 
growth time of the instability is of the order of the time it takes 
an Alfven wave to travel around the star on a field line. This 
and other instabilities were reviewed by Spruit (1999). 

In Braithwaite (2005) we investigated the development of 
this instability numerically, the results confirming the conclu- 
sions from the previous analytical work. We showed that a 
toroidal field of strength B = Bqw/wq (where Bo and wq 
are constants) in a stably stratified atmosphere is unstable on 
the axis to perturbations of the m = 1 type. We confirmed that 
the growth rate a of the instability is approximately equal to 
the local Alfven frequency loa, given by: 
B B 



where va is the Alfven speed and p is the density. We also con- 
firmed that rotation about an axis parallel to the magnetic axis 
can suppress the instability if £1 > loa- However, this stabi- 
lization only takes place for the case with neither thermal nor 
magnetic diffusion (k = r\ = 0). It is not entirely certain what 
effect rotation may have when these two are present, although 
it seems very likely that in the limit £1 ^> oja, the growth rate 
is merely reduced by a factor wa/O, so that: 

(fi«w A ), (2) 



LOA 

n 



In the unstratified case where the magnetic diffusivity is 
zero (rj = 0), all vertical wavelengths are unstable. However, 
stratification stabilizes the longest vertical length scales. This 
is because it discourages any vertical motion, which is greatest 
in modes of large vertical scale. Magnetic diffusion stifles the 
shortest wavelengths, since fluctuations in the magnetic field 
produced by the instability are smoothed out by the diffusion 
at a rate which depends on the length scale of the fluctuations. 
If n is the vertical wavenumber of the unstable mode, then 



a 2 
— > n 



> 



N 2 



(4) 



where tjjq is some measure of the extent of the field in the hor- 
izontal direction. 

The two effects can conspire to kill the instability com- 
pletely when the upper and lower limits on wavelength meet 
each other. This puts a lower limit on the field strength for in- 
stability, expressed by the inequalities: 
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(Q < w A ), 
(fi > w A ). 



(5) 
(6) 



For the core of the present Sun, this yields a minimum field 
strength of the order 3 x 10 3 G. 

1.2. Expected properties of the dynamo 

In this section we summarize the scenario of Spruit (2002) for 
the generation of a self-sustained magnetic by differential rota- 
tion in a stably stratification. 

A weak initial field is wound up by differential rotation. 
After only a few differential turns the field is predominantly 
toroidal (B^ ^ B m and B^ 3> B z ). Eventually the field 
strength at which instability sets in will be reached (given by 
Eqs. (|5} and l|6}). There are various time-scales of relevance 
here, the shortest being the reciprocal of the Brunt-Vaisala 
(buoyancy) frequency N, given by A^ 2 = (g/T)(dT/dz + 
g/cp). Provided the star is rotating at less than the break-up 
rate, the rotational frequency £1 will be smaller than N. In most 
stars, £1 will however be greater than the magnetic frequency 
wa (see Eq. Q). The differential rotation time-scale is given 
by 



(i) T dr = (tu Q d z ny 



(7) 



It is this time-scale which will determine how quickly the 
initial field is wound up into a predominantly toroidal field. 

Spruit (2002) derives properties of the dynamo in the case 
where: 



(8) 



($1 » loa)- 



(3) 



This is the most realistic regime. In addition, we expect in a 
real star to have ra r of the same order as, but in most regions 
probably greater than, fl^ 1 . 

At the time when the instability sets in, the growth time of 
the instability is so long (from Eqs. and (0) that the field is 
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still being wound up faster than it is able to decay. However, as 
the field grows, and the instability growth rate a rises, a point 
will be reached where the field is decaying and being wound up 
at the same rate - we call this 'saturation' . The time-scale r a on 
which the field component B z is wound up into an azimuthal 
component of comparable strength to the existing azimuthal 
field is given by 



(9) 



Thus, the greater B z , the shorter the amplification time-scale. 
At this point, we need to know the value of B z . This is provided 
by the instability which produces a vertical component from 
the azimuthal field by the unstable fluid displacements. If these 
have a vertical length scale / and horizontal scale zuq, we have 



B^I/zuq, 



(10) 



So, B z is greatest for the unstable modes with greatest vertical 
wavelength, i.e. with lowest wavenumber n = 1/7. If we equate 
this minimum amplification time-scale to the instability time- 
scale cr -1 , we get, using Eqs. (|4) and < l 1 01 . 



L0 A 



Nr dr 



(11) 



in the case where ^> uja- In the slowly rotating (£1 -C u>a) 
case, we find that uja drops out of the equation leaving N = 
t ( J^ 1 3. In this case, therefore, the field will either grow until it 
is strong enough to kill the differential rotation itself, or it will 



decay until the ft 3> u>a regime is reached. If r a < 
field will grow, i.e. if 

N < rZ 1 . 



the 



In real star, of course, this will not be the case except in or 
near a convective zone. Therefore, the field will decay until the 
regime ft ^ uja is reached. 

2. The numerical model 

We take advantage of the fact that the Tayler instability is lo- 
calised on the meridional plane, and model just a small sec- 
tion of the star on the rotation axis. This is similar to the ar- 
rangement used in Braithwaite (2005), where we looked at the 
Tayler instability in the absence of differential rotation, mod- 
elling only a small volume on the magnetic axis of symmetry. 

Inside the computational box, the plasma is rotating in the 
horizontal plane about an axis passing through the centre. The 
rotation speed ft is independent of distance from the axis w 
and a function of just height z. The bottom face of the box lies 
in the plane z = 0. The computational box has a height L and 
width 2L. 

We use a system of Cartesian coordinates. At first glance 
one might instead consider using cylindrical coordinates, as 
these appear to be more suited to the task. This is indeed 
the case if one is wanting to handle the problem analyti- 
cally. However, cylindrical coordinates not only make numeri- 
cal modelling more time-consuming per grid box, but also in- 
troduce special points (the axis). This coordinate singularity 




Fig. 1. The computational box. The cylinder represents the vol- 
ume where the gas is being acted on by the rotational force in 
Eq. GU- 



is a known problem in all grid-based codes in cylindrical and 
spherical coordinates. Since the phenomenon we wish to inves- 
tigate lies on the z-axis, it is better to use Cartesians so that we 
can be sure that our results are not merely an artefact of the 
code. The disadvantage is that some space in the corners of the 
computational box is wasted. When the output of the code is 
analysed, a conversion into cylindrical coordinates is first per- 
formed. The computational setup is illustrated in Fig. [2 



(12) 2.1. Implementation of the differential rotation 



We want the rotation rate of the gas to have a time-average de- 
pendence on height of the form il = Slo + zdfi/dz. To achieve 
this we apply a force F per unit mass of the form: 

F(vj,z) = (vo — v)/rf where vo = — — zvj&a, (13) 

dz 

where v is the velocity field, T{ is a time-scale, which can be 
chosen, and is the azimuthal unit vector. No force was ap- 
plied in the vertical direction. The gas at z — is therefore 
not rotating. We could in principle add a constant O to the 
rotation speed, but this would result in the gas moving more 
quickly, and the time step of the code would go down. It is bet- 
ter to include ilo indirectly: we transform to the rotating frame 
and add the Coriolis force 2v x Oo to the momentum equa- 
tion. The centrifugal force can be ignored since its only effect 
would be a change in the equilibrium state. 

This force F is applied to the gas out to a radius L, i.e. to the 
sides of the computational box. This means that the corners es- 
cape this force. This was found to be the best way to reduce the 
effect of the geometry of the box on the physical phenomenon 
of interest - the gas in the corners is roughly stationary and its 
effect on the rotating gas in the middle is minimal. 

In a star, differential rotation is driven by two mechanisms: 
magnetic braking, which acts on the surface of the star, and 
evolution, when radial shells are contracting and expanding. 
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The former would be difficult to model, as exerting a force at 
the boundaries could be expected to cause problems at those 
boundaries. Exerting a torque by means of a distributed force, 
as is done here, is a good way to model differential rotation 
caused by evolution. 



gas to be rotating as fast as possible (to maximize the chances 
of creating a dynamo) but still comfortably below the sound 
speed. With a magnetic field this weak, there is still plenty 
of space to fit the rotation time scale between the sound and 
Alfven time scales. 



2.2. Boundary conditions 

In all three directions, we employ mirror-like boundary con- 
ditions. This means that unknown quantities just outside the 
boundaries are copied from the same quantities just inside the 
same boundary. This is somewhat more difficult to implement 
than the periodic conditions used in Braithwaite (2005), but is 
necessary for the following reason. The gas in this case is ro- 
tating in a horizontal plane; with periodic boundary conditions 
the gas at the boundary would be moving in the opposite di- 
rection to the gas immediately adjacent on the other side of the 
boundary. This would cause turbulence - we are only interested 
in one kind of instability and introducing another would surely 
confuse the issue. In fact, periodic boundary conditions were 
tried at first, and it was found that a dynamo process was pro- 
duced even in the case of uniform rotation (dil /dz = 0), as 
a result of the hydrodynamic turbulence created at the bound- 
aries. This argument does not of course hold at the boundaries 
in the vertical direction, but there is a good reason for having 
the same conditions there. Sound and internal gravity waves 
can propagate upwards, their amplitude growing as they do so. 
With periodic conditions at the vertical boundaries, a wave can 
go through the top of the computational box and come back at 
the bottom, to continue its upwards travel. Its amplitude rises 
exponentially on a timescale roughly equal to A^ 1 and will 
eventually get out of control. In Braithwaite (2005) we pre- 
vented this by having gravity point in opposite directions in 
the top and bottom halves of the volume modelled, the draw- 
back being the computational expense of having to model the 
same thing twice. The method used in this study is, as far the 
the end result in concerned, the same as that previous method, 
but computational cost has been exchanged for programming 
complexity. 

2.3. Initial conditions 

We begin with uniform temperature; pressure decreases expo- 
nentially with increasing z, in hydrostatic equilibrium. The ini- 
tial velocity field we create by running the code with no mag- 
netic field until a steady state has been reached. The initial mag- 
netic field should in principle be unimportant, except that it 
must have a vertical component and must be weak. We choose 
therefore the simplest field imaginable: a uniform vertical field 
B = Boe z , where e z is the vertical unit vector. The field en- 
ergy must be weak compared to both the thermal energy and 
the kinetic energy, or in other words, the Alfven speed must be 
much less than both the sound speed and the rotation speed. Bq 
is chosen so that the ratio of thermal to magnetic energy den- 
sities j3 = 10 5 , or so that the ratio of sound to Alfven speeds 
is 240. The ratio of rotation speed to Alfven speed depends of 
course on the values we choose for fl and dSl /dz. We want the 



2.4. Free parameters 

The goal is to produce a self-sustaining magnetic field, but 
there are few clues as to precisely what conditions may be nec- 
essary. We have a fair number of parameters to play with. 

The values of flo and dfi/dz will have an effect. It is ex- 
pected that a value of Qq above oja will slow the growth rate of 
the Tayler instability. This would then increase the saturation 
field strength. Whether this makes a dynamo any more or less 
likely to appear in our model is not clear. What is certain, how- 
ever, is that a large value of dCl/dz will be conducive to the 
appearance of a self-sustaining field. We set this therefore in 
all runs to a high value but such that the flow speed is still com- 
fortably below the sound speed. We set dQ/dz — c a /5L 2 , so 
that the gas is moving at a maximum of one fifth of the sound 
speed. 

Another choice to be made is the relaxation time-scale Tf of 
the driving force. Setting it too low would inhibit the instability, 
as it would hold the plasma to too stiff a velocity field. Too high 
a value, on the other hand, may mean that the driving force is 
insufficient to make the plasma rotate in the required manner. 
Various values are used in the results reported below. 

The code contains an artificial diffusion scheme designed 
to maintain stability. It includes terms for all three diffusivities 
(kinetic, thermal and magnetic). The adjustable coefficients in 
this scheme were set to the experimentally determined mini- 
mum value needed for numerical stability. In the simulations 
presented in Braithwaite (2005), it was possible to turn off 
this scheme completely, since we were dealing with a body of 
plasma which was stationary at t = Q and whose movement 
we only wished to follow while the velocities were small. This 
is unfortunately not the case here, as the plasma is necessarily 
moving fairly quickly. 

In an ideal simulation, all unstable wavenumbers would be 
modelled between the two limits in Eq. ©. The wavenumbers 
accessible numerically are given by the vertical size of the com- 
putational box and the spatial resolution (the Nyquist spatial 
frequency): 



To see an instability at all, we need this numerical range to 
overlap with 1 141 

3. The numerical code 

We use a three-dimensional MHD code developed by Nordlund 
& Galsgaard (1995), with extensive modifications, chiefly the 
mirror-like boundary conditions described in Sect. 12.21 The 
code uses a staggered mesh, so that variables are defined at 
different points in the grid-box. For example, p is defined in 
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the centre of each box, but u x in the centre of the face per- 
pendicular to the x-axis, so that the value of x is lower by 
idx. Interpolations and spatial derivatives are calculated to 
fifth and sixth order respectively. The third order predictor- 
corrector time-stepping procedure of Hyman (1979) is used. 

The high order of the discretization is a bit more expen- 
sive per grid point and time step, but the code can be run with 
fewer grid points than lower order schemes, for the same ac- 
curacy. Because of the steep dependence of computing cost on 
grid spacing (4th power for explicit 3D) this results in greater 
computing economy. 

For stability, high-order diffusive terms are employed. 
Explicit use is made of highly localised diffusivities, while re- 
taining the original form of the partial differential equations. 

4. Results 

We present results for a number of different setups. First, we 
look at the unstratified case, both with and without rotation 
(flo = and ilo ^ 0). Then, the stratified case, again both 
with and without rotation. 

We ran the code with 64 grid points in the horizontal direc- 
tions and 32 in the vertical. This is a trade-off between the res- 
olution needed to obtain dynamo action and the large number 
of time steps needed. The time step is set by the sound crossing 
time, but the (Alfvenic) time-scales of interest to us are neces- 
sarily much longer. As explained, the high order of the spatial 
discretization of the code guarantees a high effective resolution 
even at this limited number of grid points, and self-sustained 
magnetic fields were readily obtained at this resolution in all 
cases studied. 

4.1. An unstratified case 

In the absence of gravity, there is no maximum to the verti- 
cal length scales that are unstable (see Eq. ©). Since it is the 
maximum unstable wavelength which determines how quickly 
the field is wound up, the dynamo never becomes saturated; the 
field simply continues to grow until it is strong enough to kill 
the differential rotation. 

However, the fact that we are conducting this simulation 
inside a box of finite dimensions changes things somewhat, by 
imposing an artificial maximum wavelength. This enables the 
field to find a saturation level, and the instability operates at the 
largest length scale that fits into the numerical box. 

With the unstratified setup, the production of a statistically 
steady, self-sustained magnetic field was observed. To under- 
stand the properties of the field produced, it is useful to see the 
evolution of the mean magnetic energy density, split up into 
its poloidal and toroidal components. To this end, B^/8ir and 
B 2 v j%-K = [B% + -B^)/8tt are plotted in Fig.0 The time on the 
horizontal axis of this graph is expressed in units of the sound- 
crossing time r s = L/c s . The field, initially poloidal, becomes 
mainly toroidal as it is wound up by the differential rotation. 
This happens very quickly, over the time-scale Td r - Both com- 
ponents then grow, more slowly, until the saturation level is 
reached, when the field is being destroyed by the instability at 




1000 2000 3000 

time / r s 



Fig. 2. Amplitude of the self-sustained magnetic field in an un- 
stratified case with maximum shear velocity 1/5 of the sound 
speed. Averages of B\ /8tt (solid line) and B^/Sti (dotted line), 
in units of the thermal energy density. The field is predomi- 
nantly toroidal. Time is in units of r s , the sound crossing time 
of the box. 

the same rate at which it is being amplified by the differential 
rotation. 

Fig.|3]shows contour plots of the vertical component of the 
magnetic field B z and of the azimuthal component B$, aver- 
aged in the azimuthal direction, as a function of w and z. The 
nine panels are taken at nine different times: t = 513, 664, 694, 
724, 754, 785, 815, 845 and 875r s (in units of the sound cross- 
ing time across the box). At the time of the first frame, the az- 
imuthal mean of B z is positive everywhere, as it is at the begin- 
ning of the run. B^ has been produced from the winding-up of 
this positive B z by differential rotation of positive df2/dz, so is 
also positive almost everywhere. Then the instability produces 
a new vertical component which points predominantly down- 
wards: B z switches from positive to negative. The azimuthal 
component does likewise. 

Looking at Fig.|2] it can be seen that the mean field strength 
is lower than usual during this reversal period. This reversal in 
the prevailing direction of the field happens three times during 
this run (also at t = 2200 and 3000r s ), and there is no reason 
not to presume that it should continue to happen were the run 
continued. 

A change of the field direction throughout the entire cylin- 
drical volume is not surprising when one takes into account the 
fact that only the longest wavelengths are unstable. Most of the 
unstable range of wavelengths lies outside of the range which 
can be seen in this simulation - we cannot see unstable wave- 
lengths longer than the size of the computational box. In reality, 
we would see instability over a range of length scales; in this 
unstratified case, up to infinity. 

4.1 .1 . Torques 

The torques produced by the dynamo process vary through the 
box. As an average measure of the magnetic torque, we inte- 
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Fig. 3. Azimuthal average of B z (positive values shaded light, negative dark), as a function of w (horizontal axis) and z (vertical 
axis). Contours show B^ (black lines, with ticks pointing towards negative values). The thick lines show where B z and B^ are 
equal to zero. The nine panels are taken at nine different times: t — 513, 664, 694, 724, 754, 785, 815, 845 and 875t s , arranged 
top-left, top-middle, top-right, middle-left, etc. The predominant sign of polarities of both B z and B^ changes over this period. 
Unstratified case (N = 0). 



grate the azimuthal component of the Lorentz shear stress, mul- following way: 
tiplied by the lever arm vu, over a horizontal plane. We want to 
check that the torque acting on the plasma really is magnetic 

in origin, and not kinetic, since the velocity field also produces 27r L 

a shear stress. We calculate the two respective torques in the T m (z,t) = zud(f> dza 

J 4>=0 Jm=Q 
r 2iv 



Bd>B z 
i — - — . 
4tt 



T v (z,t) = / mdcf) dzu vopv$v z . 

J ib=0 Jtu=Q 



(15) 
(16) 
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Fig. 4. Mean torque from the magnetic field T m (solid line) and 
the velocity field T v (dotted line), in units of (2/3)irPL 3 . Non- 
stratified (N = 0), non-rotating (f2 = 0) case. 

The averages over z of these two torques are plotted in Fig.|4] 
This confirms that the torque is chiefly magnetic. The torque 
from the velocity field is almost ten time smaller, sometimes 
even negative, i.e. it has the effect of increasing the differential 
rotation. 

The torques T m and T v are plotted in units of (2/3)7rPL 3 ; 
this is the torque which would exist if the shear stress were 
equal to the gas pressure P. In accretion disk parlance, this 
corresponds to a viscosity paramater a = 1. In fact, it is of 
some interest to compare the shear stress in this model to that 
found in an accretion disc, since they are caused by similar 
processes. That the stress in an accretion disc is of the order 
of the gas pressure should not be surprising, since the kinetic 
energy of the orbital motion, which drives the dynamo, is of the 
same order as the thermal energy. In a differentially rotating 
star, however, the energy available to power the dynamo (the 
difference in kinetic energy compared with a uniformly rotating 
star with the same angular momentum), is much less than the 
thermal energy. We should therefore also expect the magnetic 
energy to be much less than the thermal energy, as it cannot be 
greater than the energy of differential rotation. 

We have seen that the torque exerted by the magnetic field 
is much less than the torque which would exist if the shear 
stress were equal to the gas pressure, as we expected. A more 
interesting quantity to compare it to, in this context, is the 
torque which would exist if the shear stress were equal to the 
magnetic pressure B 2 /8ir. To do this, we calculate a dimen- 
sionless efficiency coefficient e, defined such that: 



T m (z,t) = e(z,t) J tzj: 

j>=0 



dzu i 



B 2 



(17) 



vj=0 



This efficiency e, or rather, the average of it over all z, is plotted 
is Fig. |5] After the initial winding-up phase, its value tends to 
stay at around 0.08, except for the field-reversal phases when 
the field is weaker. 

If the vertical and azimuthal component of the field were 
everywhere equal, and if the component in the m direction 
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Fig. 5. Efficiency coefficient e, as defined in Eq. (II 71 for the run 
shown in Fig.|4] 



were zero, we would have e = 1. As in the case of MRI 
turbulence in accretion disks, however, the field is mainly az- 
imuthal. Comparing Eqs. d!5l > and d!7l >. and assuming that the 
ratio £>0 / B z is the same everywhere and that £? ro = 0: 



B+Bz/An 2{B 4> /B Z ) 



B 2 /8tt 1 + {B^/B z f 



(18) 



Looking at Fig. [5] we can estimate that B,p/B z « 1.8, and 
the above equation gives us e w 0.8. The main reason for the 
low torque efficiency may then be that the ratio B^ jB z is not 
constant, rather, it is low at small m and high at large m. This 
is confirmed by looking at the last frame of Fig. [5] for instance: 
B z is strongest very close to the axis, as the Tayler instability 
is strongest there, and B^ is strongest somewhat further from 
the axis, because the winding-up effect is stronger at larger w. 

4.1.2. Parameter dependence 

For the run discussed above the net rotation Hq = 0, (the rota- 
tion rate at z = 0, the middle of the box), and the damping time 
of the applied force Tf was set equal to the sound crossing time 
r B . For values of Tf much different from this, a self-sustaining 
field was not produced. If too low a value is used, the veloc- 
ity field is too 'stiff and the instability is unable to take hold, 
although the field does reach the required strength. If Tf is too 
high, on the other hand, the differential rotation is slowed down 
by the Lorentz forces from the magnetic field before the neces- 
sary field strength for instability has been reached. 

The code was run with a range of values of T[/t s : 0.1, 1.0, 
10 and 100. Fig- EJ shows the evolution of the mean magnetic 
energy density B 2 /8ir. It can be seen that a self-sustaining field 
was observed only around Tf /r s = 1. It was also found that the 
dynamo is not produced if a lower spatial resolution is used. It 
therefore seems that where the dynamo was produced, the con- 
ditions were only just sufficient. So it should not be surprising 
that it also only works in a narrow range of the parameter Tf . 
A higher spatial resolution will be needed to obtain dynamo 
action in a less restricted range of parameter values. 
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Fig. 6. Mean magnetic energy density against time, for runs 
with N = and t { /t s = 0.1 (solid line), 1.0 (dotted), 10 
(dashed) and 100 (dot-dashed). Energy density is in units of the 
thermal energy density, time in sound crossing times r s . Only 
in a narrow range of T{ /t s is a statistically steady magnetic field 
maintained. 

4.1.3. Rotation 

In the above we have shown a working dynamo in simulations 
with differential rotation, but with no rotation overall. This is of 
course unlikely to be the situation inside a real star, so we shall 
now model a plasma with a net rotation by adding a Coriolis 
force to the momentum equation. We find that rotation above a 
certain speed is able to stop the dynamo from working, presum- 
ably because with the given resolution it is only just possible 
to see the instability working, and small changes can push the 
dynamo action just outside the effective range. We shall see in 
the next section that in the stratified case the dynamo is more 
vigorous; active in a less restricted range of parameters. 

We run the code as above, but with values of £7o which will 
put us into the £1 ^S> oja regime: Ho/t^ 1 = 0.1, 0.3 and 1.0. 
Fig-0shows the magnetic energy density as a function of time, 
for these runs, as well as for the fio = run. It can be seen 
that when the rotation speed is above some certain threshold, 
the field, although still wound up at first, decays again without 
reaching a steady state. This behavior is not fully understood. It 
is expected that rotation will slow the growth of the instability, 
and will therefore increase the minimum unstable wavelength 
(Eq. (|4j), which depends on diffusivity. This could mean that 
the field never becomes unstable. 

4.2. The stratified case 

In the stratified case, the additional parameter is the buoyancy 
frequency N. The main result will turn out to be that a self- 
sustaining magnetic field is much easier to create in the strat- 
ified case. Whereas in the unstratified case, a self-sustained 
field was only produced within a narrow range of parameters, 
with stratification the range of parameters is now much wider. 
Rotation does not destroy the dynamo, even for rotation fre- 
quency as large as the buoyancy frequency. 



) in a stably stratified star 
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Fig. 7. Mean magnetic energy density against time, for N = 
runs with forcing parameter fio/T" 1 = (solid line), 0.1 
(dotted), 0.3 (dashed) and 1.0 (dot-dashed). Energy density in 
units of the mean thermal energy density. Sustained magnetic 
fields are obtained only at the lower rotation rates. 
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Fig. 8. The means of B|/8tt (solid line) and Bp/8vr (dotted 
line), in units of the thermal energy density, for a stratified case 
N = TjT 1 . The field is even more predominantly toroidal than 
in the unstratified case (compare Fig.|2}. 

There are one or two other differences. The field builds 
up to saturation much more slowly - in the unstratified case, 
saturation was reached over one or two Alfven crossing times 
(Alfven crossing times at the initial field strength). In the strat- 
ified case, it takes at least five times longer. Also, the field en- 
ergy is much steadier; the fluctuations do not seem to be as 
large. 

Fig- ID shows the toroidal and poloidal components of the 
magnetic energy, for a run with N = t^ 1 and Tf — 10r s . The 
energy in the toroidal component is around 30 times larger than 
that in the poloidal component, a larger difference than in the 
absence of stratification (cf. Fig.|2j. 

As in Sect. 14.1.11 we calculate an efficiency coefficient e 
(see Eq. iV7\). This is plotted in Fig.|9]for this stratified run. The 
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Fig. 9. Efficiency coefficient e, as defined in Eq. illi . Stratified 
(N = t s -1 ), non-rotating case. 

value settles a little lower than in the unstratified case, at around 
0.05. In this case, however, the ratio B$jB z is much greater 
(comparing Figs. [2] and |8j, so we expect a lower efficiency e. 
The reason that the same value of e is measured despite a higher 
average of the ratio j B z must have something to do with the 
respective distributions as a function of w of the vertical and 
azimuthal components of the field. 

We can try using other values of the driving-force time- 
scale Tf , to see how robust the dynamo is. For these runs, we 
used a somewhat stronger initial magnetic field, with /3 = 10 3 , 
as opposed to p = 10 5 used in the previous runs. This helps to 
speed up the initial evolution a little, but has no effect on the 
final steady state. 

We use values of n /t s of 1, 10 and 100. In this last case, the 
time-scale Tf is much greater than other relevant time-scales, 
which will also be the case in a real star. The magnetic energy 
in these runs is plotted in Fig. ^3 A self-sustaining field ap- 
pears in all three cases, although the saturation field strength 
does depend to some extent on the value of Tf. This in fact 
does reflect reality - the faster the differential rotation is being 
driven, the stronger we expect the excited field to be. 

We have now looked at three cases: non-rotating unstrat- 
ified; rotating unstratified; and non-rotating stratified. The 
fourth combination - rotating stratified (O 3> lua and N 3> 
wa) - will exist in almost the entire non-convective zone of al- 
most all stars, and is therefore the most interesting to us. To be 
sure that we are in this regime, we shall set both N and f2n to 
the reciprocal of the sound-crossing time, i.e. to r s . 

It is found that the self-sustaining field is still produced, un- 
like when rotation was added to the unstratified case. The field 
produced is of comparable strength to that produced when ro- 
tation is absent, but oscillates in a seemingly regular fashion. 
The mean values of B^ and B z go from positive to negative, 
and back again. We saw something like this in Sect. 14.11 but 
there, the reversal of the field was a more chaotic process with, 
one presumes, no regular period. Fig.^Jis the rotating equiva- 
lent of Fig. El| Except at the highest value of Tf, the field energy 
can be seen to jump up and down in a regular way. 
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Fig. 10. Mean magnetic energy density against time, for the 
stratified case TV = r s _1 . rf/r s = 1 (solid line), 10 (dotted) 
and 100 (dashed). A self-exciting field is produced in all three 
cases. 
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Fig. 11. Mean magnetic energy density against time, for a rotat- 
ing, stratified case. Tf/r B = 1 (solid line), 10 (dotted) and 100 
(dashed). A self-exciting field is produced in all three cases, 
and a oscillatory behavior is seen in one case. 



5. Discussion and Conclusions 

We have demonstrated by direct numerical simulation a dy- 
namo process as envisaged by analytical arguments in Spruit 
(2002). It feeds off differential rotation only, and not from any 
imposed small-scale velocity field. The toroidal (azimuthal) 
component of the field is produced from the winding-up of the 
existing poloidal (meridional) component. This toroidal field is 
unstable, and produces, as a result of its decay, a new poloidal 
component, which can then be wound up itself - in this way, 
the 'dynamo loop' is closed. A weak seed field is amplified in 
this way until a certain saturation level is reached, at which the 
field is being wound up by differential rotation at the same rate 
as it is decaying through its inherent magnetohydrodynamic in- 
stability. 
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This dynamo is therefore different from the traditional Schou, J., Antia, H.M., Basu, S. et al., 1998, ApJ, 505, 390. 
convection-powered dynamo models for stars with convective Spruit, H. C, 1999, A&A, 349, 189. 
envelopes, in which the dynamo loop is closed by poloidal Spruit, H. C, 2002, A&A, 381, 923. 
components induced by an additional, non-magnetic process. Tayler, R.J., 1973, MNRAS, 161, 365. 
The dynamo found here is, in this sense, similar to the MRI 
turbulence in accretion disk l |Hawley et al. 1 996 1 which is pow- 
ered by the gradient in orbital rotation in the disk. 

A dominant factor in the case of differential rotation in 
a star is the effect of the stable stratification. It strongly re- 
stricts the types of instability that can occur in the azimuthal 
field produced by winding up in the differential rotation. In 
Spruit (1999) we have shown that the first instability to set in as 
the magnetic field strength increases by this process is Tayler- 
instability (a pinch-type instability) rather than magnetic buoy- 
ancy instabilities. 

In the simulations presented here we have considered sep- 
arately the effects of rapid net rotation and a stabilizing strat- 
ification. Perhaps somewhat surprisingly, dynamo action was 
found to set in more readily (at lower spatial resolution) in the 
stratified cases than in the unstratified case. This is the oppo- 
site of what would be expected if buoyancy instabilities were 
the dominant mechanism closing the dynamo loop. This con- 
firms the conclusion in Spruit (2002) that differential rotation 
in a stable stratification leads to a self-sustained magnetic field 
in which Tayler instability of the azimuthal field is the domi- 
nant process closing the dynamo cycle. 

The results presented were obtained only at a resolution 
close to the minimum required to obtain dynamo action. As 
a result, the field configurations obtained show only mini- 
mal structure in the vertical direction. Effectively, they cover 
only vertical length comparable with the characteristic vertical 
length scale of the process (which is governed by the strength 
of the stratification, cf. Spruit 2002). 

For further progress, simulations at higher resolution will 
allow more detailed comparison with the analytical estimates. 
An important step forward, however, would be the develop- 
ment of code more suited to the nearly-incompressible, highly 
stratified conditions relevant for stellar interiors. With the ex- 
isting codes, the time step limitations due to sound speed and 
the buoyancy frequency limit the degree of realism that can be 
achieved. 
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